---
title: "Temperature and Humidity Data"
output: html_notebook
---

```{r}
library(tidyverse)
library(sf)
library(rgdal)
library(raster)
select <- dplyr::select

# Get Albers Equal Area Conic Projection
#https://spatialreference.org/ref/esri/north-america-albers-equal-area-conic/
aea <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
```

## Combine Data

```{r, eval = FALSE}
get_cells = function(myfile){
  print(myfile)
  
  raster(myfile) %>% 
    crop(y = shapes) %>%
    # Convert to spatial polygons data.frame
    as('SpatialPolygonsDataFrame') %>%
    # Convert to Sf format
    st_as_sf() %>%
    rename(value = 1) %>%
    mutate(date = myfile %>% str_extract("2020[0-9]{4}"),
           type = myfile %>% str_extract("t.mean")) %>%
    return()
}

# Get files
folder <- data.frame(file = dir("raw_data/temperature/prism", full.names = TRUE)) %>%
  filter(str_detect(file, "[.]asc") & 
           str_detect(file, "[.]asc.aux.xml", negate = TRUE))         

# Get appropriate projection
myproj <- raster::projection(raster(folder$file[1]))

remove(folder)

#####################
# Louisiana
#####################

# Import ordinary time-invariant shapefiles for county subdivisions
shapes <- read_rds("shapes/csub.rds") %>%
  filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/hurricane_zeta_counties.rds")) %>%
  st_transform(crs = myproj)

# Get dew temperature for Louisiana
data.frame(file = dir("raw_data/temperature/prism", full.names = TRUE)) %>%
  filter(str_detect(file, "[.]asc") & 
           str_detect(file, "[.]asc.aux.xml", negate = TRUE),
         str_detect(file, "tdmean")) %>%
  split(.$file) %>%
  map(~get_cells(.$file)) %>%
  dplyr::bind_rows() %>%
    select(value, date, geometry) %>%
  st_transform(crs = aea) %>%
  saveRDS("raw_data/temperature/la_prism_tdmean.rds")

# Get temperature for Louisiana
data.frame(file = dir("raw_data/temperature/prism", full.names = TRUE)) %>%
  filter(str_detect(file, "[.]asc") & 
           str_detect(file, "[.]asc.aux.xml", negate = TRUE),
         str_detect(file, "tmean")) %>%
  split(.$file) %>%
  map(~get_cells(.$file)) %>%
  dplyr::bind_rows() %>%
    select(value, date, geometry) %>%
  st_transform(crs = aea) %>%
  saveRDS("raw_data/temperature/la_prism_tmean.rds")



#####################
# California
#####################

# Import ordinary time-invariant shapefiles for county subdivisions
shapes <- read_rds("shapes/csub.rds") %>%
  filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/glass_fire_counties.rds")) %>%
  st_transform(crs = myproj)

# Get dew temperature for Louisiana
data.frame(file = dir("raw_data/temperature/prism", full.names = TRUE)) %>%
  filter(str_detect(file, "[.]asc") & 
           str_detect(file, "[.]asc.aux.xml", negate = TRUE),
         str_detect(file, "tdmean")) %>%
  split(.$file) %>%
  map(~get_cells(.$file)) %>%
  dplyr::bind_rows() %>%
    select(value, date, geometry) %>%
  st_transform(crs = aea) %>%
  saveRDS("raw_data/temperature/ca_prism_tdmean.rds")

# Get temperature for Louisiana
data.frame(file = dir("raw_data/temperature/prism", full.names = TRUE)) %>%
  filter(str_detect(file, "[.]asc") & 
           str_detect(file, "[.]asc.aux.xml", negate = TRUE),
         str_detect(file, "tmean")) %>%
  split(.$file) %>%
  map(~get_cells(.$file)) %>%
  dplyr::bind_rows() %>%
    select(value, date, geometry) %>%
  st_transform(crs = aea) %>%
  saveRDS("raw_data/temperature/ca_prism_tmean.rds")
```

## Average to City-Date Level

```{r, eval = FALSE}

mydates <- read_rds("raw_data/temperature/la_prism_tmean.rds")$date %>% unique()


get_avg = function(ourdate){
  print(ourdate)
  
  mysubgrid <- grid %>%
    # Zoom into just that date
    filter(date == ourdate)
  
  shapes %>%
    st_join(mysubgrid) %>%
    as.data.frame() %>%
    group_by(geoid, date) %>%
    summarize(value = mean(value, na.rm = TRUE)) %>%
    ungroup() %>%
    return()
}

# Import ordinary time-invariant shapefiles for county subdivisions
shapes <- read_rds("shapes/csub.rds") %>%
  filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/hurricane_zeta_counties.rds")) %>%
  st_transform(crs = aea) %>%
  select(geoid, geometry)

# Temperature - LA
grid <- read_rds("raw_data/temperature/la_prism_tmean.rds")

data.frame(date = mydates) %>%
  split(.$date) %>%
  map(~get_avg(.$date)) %>%
  dplyr::bind_rows() %>%
  saveRDS("raw_data/temperature/la_city_tmean.rds")

# Dew Temperature - LA
grid <- read_rds("raw_data/temperature/la_prism_tdmean.rds")

data.frame(date = mydates) %>%
  split(.$date) %>%
  map(~get_avg(.$date)) %>%
  dplyr::bind_rows() %>%
  saveRDS("raw_data/temperature/la_city_tdmean.rds")



# Import ordinary time-invariant shapefiles for county subdivisions
shapes <- read_rds("shapes/csub.rds") %>%
  filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/glass_fire_counties.rds")) %>%
  st_transform(crs = aea) %>%
  select(geoid, geometry)

# Temperature - CA
grid <- read_rds("raw_data/temperature/ca_prism_tmean.rds")

data.frame(date = mydates) %>%
  split(.$date) %>%
  map(~get_avg(.$date)) %>%
  dplyr::bind_rows() %>%
  saveRDS("raw_data/temperature/ca_city_tmean.rds")



# Humidity - CA
grid <- read_rds("raw_data/temperature/ca_prism_tdmean.rds")

data.frame(date = mydates) %>%
  split(.$date) %>%
  map(~get_avg(.$date)) %>%
  dplyr::bind_rows() %>%
  saveRDS("raw_data/temperature/ca_city_tdmean.rds")

```


```{r}
bind_rows(
  read_rds("raw_data/temperature/la_city_tmean.rds") %>%
    mutate(type = "tmean"),
  read_rds("raw_data/temperature/la_city_tdmean.rds") %>%
    mutate(type = "tdmean"),
  read_rds("raw_data/temperature/ca_city_tmean.rds") %>%
    mutate(type = "tmean"),
  read_rds("raw_data/temperature/ca_city_tdmean.rds") %>%
    mutate(type = "tdmean")
) %>%
  mutate(date = lubridate::make_date(year = str_sub(date, 1,4),
                                     month = str_sub(date, 5, 6),
                                     day = str_sub(date, 7, 8))) %>%
    pivot_wider(id_cols = c(geoid, date),
              names_from = type, values_from = value) %>%
  # Calculate relative humidity by two methods
  # https://earthscience.stackexchange.com/questions/16570/how-to-calculate-relative-humidity-from-temperature-dew-point-and-pressure
  # m = 7.591386 when temperature = -20 to +50 Celcius
  # Tn = 240.7263 when temperature = -20 to +50 Celcius
  # Source from Vaisala (2013). Humidity Conversion Formulas: Calculation formulas for humidity. Helsinki, Finland.  https://web.archive.org/web/20200212215746im_/https://www.vaisala.com/en/system/files?file=documents/Humidity_Conversion_Formulas_B210973EN.pdf

  # Or use this approximation
  # Lawrence, Mark G., 2005: The relationship between relative humidity and the dewpoint           temperature in moist air: A simple conversion and applications. Bull. Amer. Meteor. Soc., 86,     225-233. doi: http;//dx.doi.org/10.1175/BAMS-86-2-225
  # 100 - 5 (T - TD)
  mutate(humidity_a = 100 - 5*(tmean - tdmean)) %>%
  mutate(humidity_b =  100 * 10^(7.591386 * 
                                   (  (tdmean / (tdmean + 240.7263))  -  
                                           (tmean / (tmean + 240.7263) ) ) ) ) %>%
  saveRDS("raw_data/temperature/weather.rds")
```



## Estimate for each Week

- I need to estimate the average temperature and humidity for TWO WEEKS PRIOR to the date of test positivity rates.

```{r}
# Our Louisiana data is measured by Wednesdays
#This sets the start day of the week to be Wednesday (7 is the default, Sunday)
options(lubridate.week.start=3)

outa <- read_rds("raw_data/temperature/weather.rds") %>%
  # Filter to Louisiana
    filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/hurricane_zeta_counties.rds")) %>%
  # Identify the week of that observation
  mutate(actual_week = lubridate::ceiling_date(date, unit = c('week'))) %>%
  # Identify the week that we want to link this observation to,
  # because two weeks later (14 days), they might lead to greater transmissibility
  mutate(week = actual_week + 14)  %>%
  group_by(geoid, week) %>%
  summarize(tmean = mean(tmean, na.rm = TRUE),
            tdmean = mean(tdmean, na.rm = TRUE),
            humidity_a = mean(humidity_a, na.rm = TRUE),
            humidity_b = mean(humidity_b, na.rm = TRUE)) %>%
  ungroup()


# Our California data is measured by Mondays
#This sets the start day of the week to be Monday (7 is the default, Sunday)
options(lubridate.week.start=1)

outb <- read_rds("raw_data/temperature/weather.rds") %>%
  # Filter to Louisiana
    filter(str_sub(geoid, 1,5) %in% read_rds("raw_data/glass_fire_counties.rds")) %>%
  # Identify the week of that observation
  mutate(actual_week = lubridate::ceiling_date(date, unit = c('week'))) %>%
  # Identify the week that we want to link this observation to,
  # because two weeks later (14 days), they might lead to greater transmissibility
  mutate(week = actual_week + 14)  %>%
  group_by(geoid, week) %>%
  summarize(tmean = mean(tmean, na.rm = TRUE),
            tdmean = mean(tdmean, na.rm = TRUE),
            humidity_a = mean(humidity_a, na.rm = TRUE),
            humidity_b = mean(humidity_b, na.rm = TRUE)) %>%
  ungroup()

# Reselt to default
options(lubridate.week.start=7)

# Save output 
bind_rows(outa, outb) %>%
  write_csv("raw_data/temperature/weather_weekly.csv")

remove(shapes, outa, outb)
```
```{r}
read_csv("raw_data/temperature/weather_weekly.csv") %>%
  # We can also calculate a Z-score, 
  # which is the average of the z-scores of 
  # our more robust humidity measure and temperature
  # to show, on average, which communities 
  # have higher temperatures and humidity,
  # since these factors impede COVID spread
  mutate(weathertrans = (scale(humidity_b) + scale(tmean)) / 2,
         weathertrans = as.numeric(weathertrans)) %>%
  write_csv("raw_data/temperature/weather_weekly.csv")
```

## Weather Weekly Visual

```{r}
out <- read_csv("raw_data/temperature/weather_weekly.csv") %>%
  mutate(state = str_sub(geoid, 1,2)) 

out %>% 
  pivot_longer(cols = c(humidity_a, humidity_b),
               names_to ="variable", values_to = "value") %>%
  ggplot(mapping = aes(x = week, y = value, color = variable)) +
  geom_jitter(alpha = 0.5) +
  geom_smooth() +
  facet_wrap(~state)
```

